{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximum a posteriori (MAP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\theta_{\\mathrm{MAP}}\n",
    "&= \\arg\\!\\max_{\\theta} \\bigg (\\prod_{i=1}^{m}p(x^{(i)}|\\theta)\\bigg) p(\\theta) \\\\\n",
    "&= \\arg\\!\\max_{\\theta}\\bigg (\\prod_{i=1}^{m} \\sum_{z^{(i)}} p(x^{(i)}, z^{(i)} | \\theta) \\bigg) p(\\theta) \\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Denote the right side of $\\arg\\!\\max$ as $A$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Derivating EM algorithm for MAP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take log to both sides and apply Jensen's inequality\n",
    "\n",
    "For convex function:\n",
    "\n",
    "$$ E[f(X)] \\ge f(E[X])$$\n",
    "\n",
    "For concave function (e.g. log in this case):\n",
    "\n",
    "$$ f(E[X]) \\ge E[f(X)]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Detour 1**: 1-D illustration for Jensen's inequality:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAADXCAYAAAA+/9QPAAAABHNCSVQICAgIfAhkiAAAABl0RVh0\nU29mdHdhcmUAZ25vbWUtc2NyZWVuc2hvdO8Dvz4AACAASURBVHic7d17WFTXuT/w7xhEYTQaFO+K\nBk0zYrxCbk24nDRRjI2RS9vUSImYiLFqRI2Y5+ivwfZoEvFCq2DSegFtTw3QaAQ0NspMWq0yGjA6\nxHiqjCIqo+KFQVFg/f7wzD5sZgZmmIEB5vt5Hp7A2rc1ZmbNu9Ze+10KIYQAERERETVLJ1dXgIiI\niKg9YzBFRERE5AAGU0REREQOYDBFRERE5AAGU0REREQOYDBFRERE5AAGU0REREQOYDBFRERE5AAG\nU0REREQOYDBFRERE5AAGU0REREQOYDBFRERE5AAGU0REREQOYDBFRERE5AAGU0REREQOYDBFRERE\n5AAGU0REREQOYDBFRERE5AAGU0REREQOYDBFRERE5AAGU0REREQOYDBFRERE5AAGU0REREQOYDBF\nRERE5AAGU0REREQOYDBFRERE5AAGU0REREQOYDBFRERE5AAGU23YnTt3sGrVKjz77LPw8fGBp6cn\n+vfvj6ioKBw+fNjV1SOiNojtBlHrUwghhKsrQebOnDmDV155BRcuXLC4XaFQIDU1FbNnz27lmhFR\nW8V2g8g1ODLVBt27dw9Tp06VGsR58+ahtLQUlZWV+Pjjj9GpUycIITB37lz84x//cHFtiagtYLtB\n5EKC2pxPP/1UABAAxJQpU8y2z507V9r+05/+1Gx7bW2tSEtLE08//bTo1q2b8Pb2FmPGjBHJycmi\nurpa2q+yslI6z+jRo0VdXZ346KOPxOOPPy48PT3FE088ITZv3iw7d0hIiHTMuXPnZNvu3bsnunXr\nJgCIESNGyLbduXNHvP/+++Lxxx8XnTt3Fv379xfvvvuuqKiokO03adIkAUB06tRJ/Otf/5Jt+/LL\nL6VrR0VF2faPSeQmWqvdEKJ5bYfJ3bt3xYcffiieeOIJ4enpKXx9fcWsWbPEtWvXLO7/97//Xfz0\npz8Vffr0ER4eHqJnz57i5ZdfFl999ZW0T0u3S0RNYTDVBpkCCgBiz549ZtuLi4vF0KFDxVtvvSV2\n7txptv3nP/+5dHzDn0mTJona2lohhBB1dXVS+aBBg0RiYqLFYzIyMqRzp6amSuV/+MMfZNfNycmR\ntv3mN7+Ryu/duyeefvppi+ceP368rKEuKysTPj4+AoB46qmnxP3794UQQhgMBtGvXz8BQPj5+bGx\nI2qgtdoNIZrXdgjxMGB76aWXLO7r5+cnbt68Kdt/165dolOnThb3VygUYvfu3UKIlm+XiJrCYKoN\n6t+/v9VeVlN27twpHTtu3Dhx9uxZodfrxTPPPCOVf/bZZ9L+pjJTb1Kr1YrLly+LKVOmSNueffZZ\naX+DwSA8PDwEABEeHi679uzZs6VjfvjhB6l8/fr1UvmSJUvEnTt3xJIlS6Sy1NRU2Xl27dolbVu5\ncqUQQoioqCgBQDzyyCPiH//4h13/JkTuoDXbDSHsbzuEEOKPf/yjtC0yMlJcu3ZNJCQkSGULFy6U\n7R8SEiJ69OghPDw8RE5OjqiurhZbtmyRBT1CtE67RNQYBlNtkKenp/SBNhqNdh37H//xH9Kx+/bt\nk8oPHDgglYeEhEjl9XtjWVlZUnlRUZFU3qtXL9k1Jk6cKAAILy8vUVVVJYR42FMdMGCAACCCgoJk\n+z///PPSuS5duiSEEKK8vFwqe/nll81ex5tvvikAiC5duoj/9//+n7Tvhx9+aNe/B5G7aM12Q4jm\ntR2hoaHStqNHjwohHo5Ge3p6il69eokXX3yx0XrW1taKO3fuyAI5k9Zol4isYTDVBnXt2lX6QFdW\nVtp1rOkWGQBhMBikcoPBIJU/9thjUnn9BrH+/nfv3pVtq2/r1q1S+d69e4UQQhw7dkwqW7dunWz/\nnj17WhxKN/3069fP7HXcvHlTDB48WLZfcHCw7FYDEf2f1mw3hGhe21H/Og1v6Vmzfft2ERQUJLp3\n726x/TBpjXaJyBo+zdcG9e3bV/q9vLzcrmNv3bol/f7YY49Z/P3OnTsWj+3du7f0e9euXa1eY9q0\naejSpQsAIDc3FwCwZ88eAMAjjzyCX/ziF7L9rV3P5Pr162ZlPXr0wH/+53/KykxPJBGROVe1G4Dt\nbUf963Tr1q3Jeq1fvx6/+tWvUFBQIF1foVBY3Lc12iUia/jN1AaNHz9e+r2oqMhs+6lTpzBy5Egs\nX74chYWFsm09e/aUfq/fWNT/vUePHg7Vr0ePHpg0aRIAICcnB8D/NVphYWHo16+fbP/u3btLv9+9\nexfi4Yio9HP//n2za9y6dQsfffSRrGzRokWoqalxqO5EHVVbbzcAeQDVVDADAH/4wx+k35OTk2E0\nGlFdXW1x39Zol4isYTDVBr3xxhvS72lpaWbb161bh+LiYvz2t7/FokWLZNvGjRsn/V4/2/E///lP\n6ff6jW5zmXp5er0eubm5OHnyJABg+vTpZvuOHDlS+v1//ud/bDr/22+/jXPnzqFLly5YtmwZgIev\n4f3333e06kQdUntoNwICAqTfdTodAODy5cvw9PTEY489hmeffVa2/8WLF6Xf58yZA29vb7NAsL6W\nbpeIrHLdHUaypqamRvYUza9//Wtx6dIlcefOHZGUlCQUCoUAIDp37myWi+nPf/6zdNyECROEXq8X\nZ8+eFWPHjpXK//rXv0r7w8rchqa2VVZWCm9vbwFADB8+XAAQXbt2Fbdu3TLbt/5TMz//+c/FtWvX\nxJ///Gfh5eUlBg0aJH72s5/J9k9LS5P2T05OFkIIMXPmTKksMzOzWf+uRB1Za7YbQjSv7fj9738v\nlb/++uvi+vXrYsGCBVLZ/PnzZfsPGTJENjFer9eL5557TpYuoX5+qpZsl4gaw2Cqjbp06ZIICAiw\nOjnykUceEdu3b7d4bHR0tNXjpk+fLtu3ucGUEOZ5aawl0qyurhZBQUEW69OlSxeRn58v7Xvy5Elp\nIu3LL78s6urqhBBCVFVViVGjRgkA4tFHHxVnzpyx6d+RyJ20VrshRPPajvv374vnnnvO4jUs5Zla\nvHix2X5PPPGEeO+996S/+/TpIzumJdoloqYwmGrD7t27J9avXy9+/OMfCx8fH+Hh4SH69+8vfvnL\nX4rjx49bPa6mpkakpKSIcePGCS8vL6FUKsXTTz8tPv30U7On4RwJpv72t7/J9snOzrZap9u3b4v3\n339fDBs2THTu3Fn07NlTvPrqq9Lj0UIIYTQahUqlEsDDR6pNjyubFBcXC6VSKQCIUaNG2f34N5E7\naI12Q4jmtx137twRiYmJUlvQp08f8dZbb4mysjKzfY1Go3j33XdF3759xaOPPioiIiLEpUuXRFlZ\nmQgLCxP9+vUzy2Xl7HaJyBZc6JiIiIjIAZyATkREROQABlNEREREDmAwRUREROQABlNEREREDmAw\nRUREROQABlNEREREDmAwRUREROQABlNEREREDmAwRUREROQABlNEREREDmAwRUREROQABlNERERE\nDmAwRUREROQABlNEREREDmAwRUREROQABlNEREREDmAwRUREROQABlNEREREDmAwRUREROQABlNE\nREREDmAwRUREROQABlNEREREDmAwRUREROQABlNEREREDmAwRUREROQABlNEREREDmAwRUREROQA\nBlNEREREDmAwRUREROQABlNEREREDmAwRUREROQABlNEREREDmAwRUREROQABlNEREREDmAwRURE\nROQABlNEbkav17u6CkREHQqDKSI3odVqMX/+fGRmZrq6KkREHYqHqytARC1Lr9cjPT0dBoMBMTEx\nCAwMdHWViIg6FI5MuSmFQiH7seTYsWOYMmUKfHx84OnpiTFjxgAALly4gMjISJuvFRUVhdLS0mbX\ng5rHaDQiIyMDSUlJUKlUSElJYSBF1AhTm9dS7RLbzo6LwVQHkZqaioEDB5p9wOr/9OnTx+y4vXv3\noqyszKz80KFDeOGFFzBo0CCUlJRAo9HA398f+fn5CAsLw5IlS2yu2+LFixEaGopvvvnGbFtZWRn2\n7t1r34ulJuXl5WHBggWorKzE6tWrERUV5eoqEbVp9ds8k+eee86sfWxuW8u2s4MT1GFUVFQIAKJv\n375SWU1NjaiqqhLff/+9CAsLk8oBCGv/+2tra4W/v7/w8vISt2/flsq//fZb4e3tLbKzs+2uW1ZW\nllAqlaKwsNDi9sbqQ7YrKSkRSUlJYunSpeL06dOurg5Ru9CwzWuqPbKnrRWCbac74MhUB1JUVAQA\nGD16tFT2yCOPwMvLC0qlEqNGjbLpPAcPHsS///1vhIaGonv37gAAIQRiY2MREBCAadOm2V23iIgI\nqFQqxMbGQghh9/HUuPq39AIDA7F69WqMHDnS1dUiahcstXmNsaetZdvpHhhMdSCmD7hpblN9gwYN\nQkpKik3nOXDgAAAgNDRUKsvNzUVRURGio6ObXb/o6GgUFhZi3759zT4HmdNoNFiwYAHKy8uxevVq\nhIeHu7pKRO2KpTavMfa0tWw73QODqQ6ksLAQgLy3BACvvfaaTcdrtVoEBQVhzZo1AIClS5dK9/+z\ns7MBAGPHjpUdc/LkSfTv3x8KhQI7duwAAOzevRsKhQIDBgyQ7Wuq1xdffGHnKyNLDAYDVq5ciczM\nTCQkJGDRokXw9fV1dbWI2g1rbV5T7Glr2Xa6B7cKpoxGo6ur0KIsDT2XlZWhtrbWpuMDAwNRUFAg\nfZAvXrwIIQTKy8tx/PhxAEBAQIDsmNGjR+PkyZPo3Lmz1Gh4enpCoVDgo48+ku1ruu1kOhc1j9Fo\nRFZWFhYsWCA9pcdbekT2s9bmNcWetpZtp5tw5YSt1pSfny+SkpJcXY0W8+DBA9GlSxdpMmL9nw8+\n+MBsf1iZtHj9+nUBQPj4+MjKH3vsMQFA3Llzx+L1p06dKry9vcW5c+fEkCFDxNatW832uXnzpsVz\nN1Yfkjt9+rSYN2+eSEpKEuXl5a6uDlG7Z6nNa6w9sretZdvpHtxmZCowMBAGgwFardbVVWkRZ86c\nQXV1NQICAiCEkH7mzJmDcePG2Xwea3MB7ty5AwDw8vKyeFxsbCyqqqrw0ksvISUlBbGxsWb7KJVK\n2bnIdkajEWvXrkVycjIiIyOxfPly3tIjcoLG5j9ZYm9by7bTPbhNMKVUKhEZGYn09HRXV6VFmO7h\nN7wv7+HhYVbWGGsNi2kegbX5BD/60Y8AACqVClOnTm30GkwwZx9TzqjevXsjJSUFISEhrq4SUYdh\nbzBlb1vLttM9uE0wBUD6EtLpdC6uifNZaxBSUlIwfPhwh8/TrVs3AMD9+/fNjjlx4gRmzJiBqVOn\n4sCBAzAYDBbPfe/ePQCw6dFjergMzMqVK6FWq5GQkICYmBiph0pEzmFvMGVvW8u20z24VTAFAJGR\nkcjKynJ1NZyusQahpqYGP/nJTxw6j5+fHwDg+vXrsvKvv/4aERERyMjIwNy5c/HgwQPs3LnT4rmv\nXbsGABg2bJhNdXFXpgnmpmVgmDOKqOU4K5gCLLe1bDvdQ5sKpioqKhAXFycru3LlCqKjo6HRaKSy\nhISEZl8jMDAQJSUl0Ov1zT5HW2TpAy6EwN27d/GnP/1J+kA3pqamBjqdDh4eHmZf3qY13c6dOyeV\n7dmzB5MnT8bGjRuhUqkQFhYGHx8fbNu2zeL5S0pKZOciczqdDsuWLYNOp+MyMEQtrLE2zxp721q2\nne6hzQRTFRUVmDhxIl555RWp7F//+hc+/fRT5Obmyt7oN2/elPJy2EupVCIkJAS5ubkO17ktWLRo\nEQYMGICrV68CAPr16yetD9WpUyd4e3sjPj4eYWFhTZ7r+++/R3V1NZ588kl06dJFts20OOfRo0el\nslmzZuH+/fvSOlQeHh6YOnUqioqKMHz4cNy6dUt2joKCAgDA66+/3vwX3EHVn2AeExPDCeZEraCx\nNq+h5ra1bDvdhOseJJSbNm2a2Lx5s1m5wWAQffr0kZWdPn1avPDCC82+Vnl5uYiLixOVlZXNPkd7\nBwuP0+7YsUMAENOnTzfbv66uTowbN06MHz++2dc0HV9XV2dTfdyFWq0WcXFxYtOmTW79niRqbdba\nPGe2R2w73YPDI1PHjh3DlClT4OPjA09PT2no88KFC1JE3pSioiL88MMPePvtt822nT592mz4dcSI\nETh16pT0d1RUFEpLSy2eu+Fq3gDg6+sLPz+/DpsmwR45OTm4fPkygMbnAigUCmzbtg1nzpxBZmam\n3dfZtWsXzp49i23btsmeSLl8+TJycnKaWfv2rWEG8zlz5nCCOVEramq+VP32sbnYdroJRyKxgwcP\nis6dO4vZs2eLW7duiSNHjohp06aJQ4cOiccff1wcOXLEpvMsWbJEbNy40eK2jRs3irlz58rKjEaj\n6Nq1q/T3kSNHhL+/v9BoNGbHl5WVib1795pF7/n5+SIxMdGm+nVEaJBsTgghXnnlFQFA7N+/3+px\n+fn5wt/fX6jVapuvdejQIav/fyzVwx1kZmaKN954Q3z++eeurgqR27LW5rVEu8S2s2Nr9v+B2tpa\n4e/vL7y8vMTt27el8m+//VZ4e3uL7Oxsm8/14x//WBw9etTitnfffdcs0Dpz5owYPHiwrCwrK0so\nlUpRWFho8TyW3nAzZ84UJSUlNtezo+vbt68AIK5cudLofnq9XkRERNh83oiICHHx4kVHq9chnD59\nWiQmJoqkpCS+94hczNY2z1nYdnZcCiFsWIjIgr///e94+eWXER4eLk3mFkJg3Lhx8PT0xLFjx2w+\nl6+vLw4fPowRI0aYbQsNDcWHH34oS1R49OhRJCYm4tChQ7J9g4KCUFNTgxMnTpglNzP9Xf/lpqen\nw2g0Ys6cOTbXlag5jEYjsrOzkZ+fj6ioKISHh7u6SkRE5CTNnjN14MABAA+DHZPc3FwUFRUhOjra\nrnNVVFRg0KBBFrdZmjN19OhRPPfcc2b7RkdHo7CwEPv27bPpuiEhIVw4klqcVqvFggULUF5ejpSU\nFAZSREQdjN3BlFarRVBQENasWQMAWLp0KRQKBfr06SOtfN0wpf7JkyfRv39/KBQKKaXB7t27oVAo\nMGDAAHTr1g0eHh6yY373u99h4MCBuHbtGsaNG4e0tDQAD0eWNm/ejF/+8pdmdTOt4P3FF1/Y9Fr8\n/PzQu3dvqNVqO/4FiGxjmmCenp6OhIQELFq0qFUmmBuNRhQXF0s/eXl5yMrK4vuc2jW9Xo/i4mJo\nNBpkZWVBq9WiuLjYatZwotbU7Nt8gwcPRmlpKS5evCiNKo0dOxZFRUW4dOkSBgwYINvfYDBg4MCB\nmDJlCrKzs5GXl4dXX30V27dvR1paGrZs2SKtUdSYtLQ0HDx4ELt27TLbduHCBfj5+WHChAlmT+pZ\nus0HPFz3TKfTYdGiRXa9fqLGmAKY4OBgREZGNiuIMhgMUuZjQL4MksFgkH2J6PV6VFVVSX97eXlh\n6NCh0t9+fn7w8/Pjun7ULuj1emg0GpSUlKC4uFgqV6lUsv2MRiMuXLgAAOjduzeGDh0KlUol/Zeo\ntXg0vYu5GzduoLS0FD4+PrLbc6Y39aOPPmp2jK+vLyZPnoz9+/fj/PnziI+Px5YtWzBjxgzU1tZi\nzZo1+Oyzzxq97t69e7F161bs37/f4vYePXoAAM6fP2/zawkODkZGRgaMRiMfSyeH6fV6pKWlQQiB\n5cuXw8/PD3q9XvpsGI1GWfZ9vV4Po9Eo/V3/iwN4+AVRP3mn6Quid+/eUKlUsm1+fn58D1O7ZTAY\noNVqodFoUF5ejsDAQKkjYssKDjqdDnq9HiUlJcjLy8O1a9egUqkQGBiIwMBAJsGlFtWsYMpabo47\nd+4AeNgrtiQ2Nha7d+/GSy+9hJSUFGmF7NjYWGkl7saYRrR69uxpcbvpi8RUD1solUpMmDABGo2G\nc1msMDVSADBy5EibGraOqH6g03Bk6LvvvsO1a9dw8+ZN9OjRAxUVFVi2bJm0fciQIdL709vbW/o3\nHDJkCIYOHQpvb29pX67DR+7EYDAgIyMDWq0WEyZMQGRkZLOWTRk5cqTss2M0GqHVaqHT6ZCZmQml\nUomgoCAEBwe7bRtGLcepwZTpVlrDJ+lMTLfxVCqVFEiZrF+/vsnrbtmyxab6Wbu+NSEhIcjKymIw\n1YDRaMTmzZtx/vx5DB8+HHfv3kVmZiZCQ0MxY8YMV1evWerfFgAernlluj1WVVUlrYEFmN9mA+S3\nGUyjRmVlZTh37hwef/xxTJ8+HT169LC5N03kzky3wydNmoTZs2c7dWTVtHSY6da2KbBKTk4G8PDp\n70mTJnHEipyiWXOm3nrrLWzbtg1bt25FbGysVO7j44OKigrcvXsXXbt2lR1z4sQJvPPOOxg0aBBy\nc3Nx6dIlp7+JKysr0b17d/Tq1cvsS9DanCmTuLg4rFixgl+A9aSmpuL+/fuIioqSRk6qqqqwYcMG\nPP300y5dhNfZ84lMr69hEOTr62v1fWowGJCWlgaDwYCYmBguQkpkI51Oh4yMDAghMGfOnFZvd/V6\nPXJzc6HVatGnTx+Eh4djwoQJvE1OzebUkSk/Pz9UVFTg+vXrGDhwoFT+9ddfIy4uDnl5eSgtLcXu\n3buxc+dOvPfeew5U3Zzpy3XYsGF2HxsSEgK1Wo2YmBin1qm90uv10Gq1SEpKkt2C8vb2xjvvvINV\nq1YhJCTE4YC4fqDT2Hyiqqoq2TbAtfOJsrKykJeXh+DgYCQkJDh07oqKCixevBh/+tOfEB8fj+zs\nbNy+fRsjRoxAcnKybPHvhIQErF271hkvgcglsrKykJWVhYiICJd1yPz8/KT8gmq1GgUFBUhLS0Nw\ncDCCg4Mt3mo3GAzQ6/UwGAyc4E5m7A6mampqoNPp4OHhYfaGCwwMRGFhIc6dOycFU3v27EF0dDSy\ns7OhUqkwYsQI+Pj4YNu2bU4Ppky3aJozQhASEoKVK1cymPpfBQUFGD16tCyQMunVqxcGDx4Mg8Eg\nBSyNzSdqOGrUcNSwPc0nqt+jNk0wd0RFRQUmTpwoPU167do1FBYWonfv3sjMzERsbCzKysqk/W/e\nvIkdO3bgzTffdOi6RK6QlpaGkpISbNiwoc3cXjPdCjQYDFCr1UhLS4NSqUR4eDiCg4MBPAy4MjMz\n0atXL/j4+CA3NxfDhg1z+q1Jar/svs136tQpPPXUUxg1ahS+++472bZ9+/YhPDwcn3zyCRYvXgwA\n6NOnDwwGA44dO4agoCAAwMyZM7F161b4+/vj+PHj0lN4jvrkk0/w/vvvY9++fZg4caJsW1O3+QAg\nMTERUVFRbnu7pv58IrVajW7duuHVV1+1uO+aNWtw7tw52fw0S/OJTOoHQe1xPlFLZTCPiIjApEmT\n8M4771i8pkqlks3x0ul0mD17Nr755hunXJ+oNRiNRmRkZKCkpATLly9v8wGIWq2GRqOBXq/Hk08+\niZqaGsyYMUM23SEzMxO3b9/G8uXLXVxbagvsDqZ27tyJN998E9OnT5cScJoIITBhwgQoFAqXZBYf\nP348FAoFtFqtTcvJNNRRck7ZM5+o4aP4pvlEBoMBEyZMwGuvvWbxGuvWrcOkSZPw/PPPO7n2bY9G\no0F6ejpUKhXi4+OlL4Jjx44hKSkJhw8fRmVlJVQqFb788kssXLgQWVlZTZ63qKgI06dPx3fffWfx\noYkjR47gb3/7Gz7++GOp7MGDB+jTpw8qKioQFRWF9evXW109oOE5m5lSjsghRqMRv/3tb6XR3LYe\nSNWn1+vx29/+FkuXLkWvXr3Mtq9YsQK/+tWv3LYDTv/H7tt81uZLAQ8b723btuH5559HZmZmq94P\n37VrF86ePYvDhw/LvkQuX76MEydO2HSOtpZzypXzifR6PTZu3GgxmKqqqsLNmzctrqXYkdSfYJ6Q\nkCAbXTt06BAmTpyImTNnoqSkBDqdDkuWLEFYWBh27txp0/l37tyJd99912IgJYTAnj17zHq9Dx48\nwL179wAAixcvRmhoKLZu3YoXX3zR7BxlZWU4ceIEpkyZYs/LJnIaUyA1ZMgQxMTEtIl21V5du3a1\nGEgBwDPPPAOdTsdgipwbTAEPl3TJyclBXFwc+vTpI91zbkn5+fn44IMPkJubi6eeekq2rWEm9sa0\nVM6p9jifyM/PD8OGDUNubi4mT54s22bK7N1W5jy0hPqTZBtOMK+rq8Pbb78NDw8PfPLJJ+jevTu6\ndu2KEydOYMeOHXj22Wdtusbhw4ctdjhMSyb9+te/Rvfu3WXbSktLpX/3Z599Fh9//DHCw8Pxz3/+\n0+wz2b9/f6u3aYlaw+bNmzFkyJB2u5i80Wi0GkgBD59gb7jaBrknpwdTwMMJfQcPHsTChQtbJZj6\n/e9/j/z8fIu3O+y9tRESEoL09HSzYMqe/EQNH8UHrM8nUqlUbXY+UUxMDNLT07F69WoMHDgQd+/e\nRWlpKUJCQlyaFqEl6XQ6pKWlwdfX1+ok2YMHD+Lf//43wsPD0b17dwghEBsbi4CAAEybNs3ma505\ncwaPPfaYrKyiogJbt27Fr371K4uNeEVFBfz9/aW/IyIisGrVKsTGxuLEiRN251gjailpaWkoLy9v\n13OKlEqlWVte340bN3D58mVoNJpW+a6jtsvuYOrKlSs27TdkyBCb5o04gz3XsWU+0c2bN5GYmCgL\nngDr+Ym8vLwQGBhoc36i9kKpVGLOnDnQ6/XQ6XTw9vbG0KFD20yw50ymLMynT59GTExMo2vYHThw\nAAAQGhoKAMjNzUVRUZFsbpMtKioqZB2AI0eOID8/H4sXL0bnzp0tHnP06FE899xzsrLo6GgsXbpU\negCEyNUyMzPbzWTzxvj5+eH69eu4fv26xc7NyZMn8cwzz0hP+0VFRTGoclPNXujY1VpiPhHwMAjS\n6XS4desWXnvtNa535gZsXZRYq9Vizpw5OHHiBOrq6qTyrl274t69e/jqq6/w8ssvy445efIkJk6c\niCtXriAjIwNvvvkmdu/ejddffx0K0Z5uiwAAIABJREFUhQLV1dVS4GR68tWkYfJZIQRGjRqFv/71\nrxg1apRUbgqi3nnnHWzevNms3rY8fEHkLKYs46tXr273HUrg4ZN9R48exfTp02XTKtLT02VP8+l0\nOmRlZcFgMDCockNtMpgyGAzQaDTNnk8EwKH5RAaDAQsWLMBnn33GQKoDa24W5sGDB6O0tBQXL17E\noEGDMHbsWBQVFeHSpUsW5+gZDAYMHDgQU6ZMkdaXfPXVVzF8+HB8+eWX0jJLTUlLS8PBgwexa9cu\nWfmFCxfg5+eHCRMmWJy/wWCKWovRaMSyZcs63IoAeXl5yM3NxYgRI+Dt7Y2zZ8+iW7duFhP2Mqhy\nT202mFKr1QBcl58oOTkZI0eO5G2TDsiU86agoMDunFE3btyQEvddv34dwP8to3Tnzh1069bN4nGv\nv/46Dhw4gFOnTiE0NBQffvghAOCf//wnPvvssyavu3fvXqxcuRL79+83W+j71q1b6Nmzp6xO9TGY\n6njqZ803uXLlCubNm4d58+bJvsBbM2v+2rVr4eXl1W4nnDfGNN3BaDSaLapsCYMqNyPIooKCAjFv\n3jxXV4OcLDc3V8TFxYlNmzaJyspKu48/ePCgACDCwsKkMg8PDwFA1NTUWD3ub3/7mwAghg0bJr74\n4gupfMGCBTZd96233hLXr1+3uO3BgwcCgOjcubPF7QAEP+odx40bN0RQUJD47//+b6nsyJEj4sMP\nPxTe3t7CYDDI9n/rrbdERkZGi9crNzdXzJs3r1mfq47s9OnTIikpScybN0+o1WpXV4daSJscmWor\n5s+fj/j4eJcvYUKO0+v1SE9Ph9FoRExMTLP/n65fvx4LFy7Ee++9h3Xr1gEAPD098eDBA9TW1qJT\np04WjysuLsbIkSMxefJk5OTkNPt1WFJTU4POnTvD09MT1dXVZts5MtWxWMuaf+3aNQQEBODq1auy\n8tbImm8wGLBs2TKnLLHUUanVamRlZcHX1xeRkZH8XulgLLf8BOBhEk/T7UZqn0y39JKSkhAYGIjV\nq1c71IhZSg1iurV3//59i8ecOHECM2bMwNSpU3HgwAHZJHNnMCXxbJiTitqeY8eOYcqUKfDx8YGn\npyfGjBmDCxcuIDIy0qbji4qK8MMPP+Dtt98223b69GmL7+0RI0bg1KlTAICoqCiUlpZaPb9CoZD9\n2CojIwPBwcEMpBoREhKClJQUBAcHIzk5GStXrjR7OIraLwZTjQgPD8c333zj9C8/ah0ajQYLFixA\neXk5Vq9e7ZT5b5aCKdMXiKX5Sl9//TUiIiKQkZGBuXPn4sGDBzZnSLeV6cGMYcOGOfW85FyHDh3C\nCy+8gEGDBqGkpAQajQaPPvoowsLCsGTJEpvO0VjW/NOnTyMgIMCs3FLWfGujVGVlZdi7d68dr+rh\nyNfp06dtDgjdnSmoUqlUSEpKklZZoHbOxbcZ27xNmzaJ7du3u7oaZIeSkhJpjsLp06eddt4HDx6I\nLl26CA8PD3Hv3j2pfNasWQKA0Gg0sv13794tPD09xd69e6XjfXx8xJgxY5xWJyGEOHTokAAg4uPj\nLW4H50y5XG1trfD39xdeXl7i9u3bQgghvv32W+Ht7S2ys7NtPs+Pf/xjcfToUYvb3n33XbFx40az\n8jNnzojBgwdLf2dlZQmlUikKCwutXsee90xiYqLIz8+3aV+Sq6ysFJ9//rl44403RGpqKuebtWMc\nmWpCVFQUNBqNlLeK2i6j0Yi0tDQkJSVBpVIhJSXFqfMSvv/+e1RXV+PJJ59Ely5dpHJTj/zo0aOy\n/WfNmoX79++jT58+AAAPDw9MnToVRUVFGD58OG7duuWUehUUFAB4+MQgtU2mrPmhoaFOz5pvYm1k\nylLWfJVKhdjYWIfn0eXl5UEI0WiSW7JOqVQiKioKGzZsgBACCxYsaLVk1+RcDKaa4OvrCz8/P2g0\nGldXhRqRl5eHBQsWQAiBlJSUFlnuxtpSShMnTsS4cePwl7/8RVZeXl4OIQSCgoKksi1btkAIgf/5\nn/9Bjx49nFKvv/zlLxg/fjxeeeUVp5yPnM9a1vzo6Gi7ztMwa3591uZMWcuaX1hYiH379tl1/fqM\nRiOysrIQExPT7HPQQ76+vpgzZw6WL18OnU6H+fPn8zunnWEwZYPIyEjk5eW5uhpkgU6nw7Jly6BW\nq5GQkIA5c+a0WKJVa8GUQqHAtm3bcObMGWRmZrbIta3ZtWsXzp49i23btpnNo7l8+bLTnxwk+2i1\nWgQFBWHNmjUAgKVLl0KhUEjB/tixY82OOXnyJPr37w+FQoEdO3YAAHbv3g2FQoG6ujp4eMhXAfvd\n736HgQMH4tq1axg3bhzS0tKkbeJ/F83+5S9/KTtm9OjRAIAvvvii2a8tLy8PEyZM4FNpTuTn54fl\ny5cjJiYGmZmZWLlypSx5NbVhrrvD2L7MmzeP8wLakPLycpGamipmzpwpcnNzW+War7zyigAg9u/f\nb3F7fn6+8Pf3b7VcMocOHRL+/v5mc7VM8L/zXsA5Uy43aNAgAUBcvHhRCCHEmDFjBABx6dIli/uX\nl5eLzp07i2nTpgkhHuZwUigUYsSIEeL777+3+bqpqakiOjrarFyv1wsAYsKECRaPa+o9U1lZKeLi\n4kR5ebnNdSH75ebmipkzZ4rU1FS3/be+ceOGmDlzphBCiNmzZwtfX1/RpUsXMWrUKLO2eOHCha6o\nohBCCLawNsrPzxdJSUmurgYJITIzM0VcXJzYvn17q07Y7Nu3rwAgrly5YnUfvV4vIiIiWqU+ERER\n0pcztV3Xr18XAISPj49U9thjjwkA4s6dO1aPmzp1qvD29hbnzp0TQ4YMEVu3bhVbt24Vs2bNsum6\nX375pXj66adFRUWF2babN2+a1am+poKpzz//XGzatMmmepBjKisrxfbt28Ubb7whMjMz3WqSesME\ntZGRkeLSpUuiurpa7Ny5U/Tv31+2f2slqLWEwZQdnP10GNnHlJU+KSlJlJSUuLo6RDbpaFnzOSrl\nGuXl5W6XSX3atGli8+bNFrdVVlbKnlIV4mG2+RdeeKE1qmaGc6bsEBkZySctXMBgMGDlypVIT09H\nTEwMsyxTu2Jprp1pfltjiTFNC2CrVCpMnTpVKl+/fr1N192yZQt8fHwa3ceexJwmprlSvr6+dh9L\nzefr64vly5cjPj5emk/VlpN+tmSCWuDh3MJf/OIXsjJbE9Q2NzltYxhM2SEkJAQGg4ETAluJKXt5\nYmIiVCoVVq1a1aFWoif30JGy5huNRuzbt4+pEFxo5MiRZkk/21rqnpZOUCuEwJ49e7B8+XJZua0J\napuTnLZJLhkPa8c4d6p1mOZFbdq0ibcTqF0bN26cACBOnDghlY0dO1YAEKWlpWb7//3vfxd+fn5C\np9OJr776SgAQ69atc2qdzp8/LwCIwMBAi9th5Tbf559/zvavDamsrBSbNm0ScXFxrfYgTlNaOkFt\nXV2dSE1NtfjZsTdBrbX3eXNwZMpOHJ1qWVqtFvPnz4dOp5NSHfB2ArVXNTU10Ol08PDwkKUQMI2w\nnjt3Trb/nj17MHnyZGzcuBEqlQphYWHw8fHBtm3bnFqvkpISWT1stW/fPi4b04YolUrMmTMHCQkJ\nUKvVUtvpSi2ZoLaiogLr1q1DdHQ0Bg4caHZMSyaobQqDqWaIjIyU5XIhx+l0OmleVGRkJJYvX878\nNdTudaSs+Wq1Gn5+fvxctkEjR47E6tWrERkZieTkZKxdu9Zl6/21VILaI0eOIC0tDfPmzUOvXr0s\nHtNSCWpt4pTxLTfEvFPOUT9f1Oeff+7q6hA51Y4dOwQAMX36dFl5XV2dGDdunBg/frxL6mW6dl1d\nncXtsHD7g21e+2BKpRAXFycyMzNb7boFBQUiMDBQdOrUSZbfrmvXrgKA+Oqrr2T7FxUViX79+gkA\nUjqDL774QgAQ/fv3Fz169BD379+X9vf19ZWdt1evXrLz1dXViZEjR4rvvvtOVp6XlycAiHfeeces\nzpbe583Fkalmio+PR1ZWVpub+NdemCaXL1iwAN7e3i22BAyRK3WUrPlarRYAOPG8HVAqldJTz6al\naVrj1l9gYCAKCgowYMAAAMDFixchhJCeSm24buTo0aNx8uRJdO7cGdnZ2QAAT09PKBQKfPTRRwgI\nCJDdBjctz2X6uXbtmux8mzdvRkBAAEaNGiUrN42kHj9+3LkvuCGnhGRuas2aNWL79u2urka7UllZ\nycnl5DY6Stb8NWvWtJkJzmSf/Px8MXPmTJGcnNzi7W1zEtRaSk4rhGiVBLUN3+eOYDDlgPLycjFz\n5kwmkLSRWq2Wkm4y+Sm5g46QNd/UzrlT5u2OprVu/TUnQa215LRCtHyCWmcGU4r/PSE1U2ZmJoqL\ni83yXdD/0Wg0yMzMhK+vLyIjIzmBlagdSU9Ph9FoxJw5c1xdFXKQXq9Heno6DAYD4uPjnd4Wr1+/\nHgsXLsR7772HdevWAXh46+7Bgweora1Fp07mM4uKi4sxcuRITJ48uUUWZq+pqUHnzp3h6emJ6upq\n2TbTbW5nhEGcM+WgqKgoGI3GVp/30B5oNBrMnz8fmZmZ0j18BlJE7YtGo8HkyZNdXQ1yAj8/Pyxf\nvrzFnvqzN0FtSyenBZqfoNZeDKacYNGiRcjLy2vTqf0dYTAYoNVqodVqbZpwb0pzkJmZicjISKSk\npDBzOVEzFBcXy35amykdApdv6lhCQkKQkpKC3r17Y9myZcjLy3PKeS0FU6b3zvXr12X7fv3114iI\niEBGRgbmzp2LBw8eYOfOnU6pR32mierDhg1z+rnr82jRs7sJX19fxMTEIDk5GatWrYJSqXR1lZwm\nPT0darUagwcPBgCsXbsWb775psWeqk6nQ1ZWFgwGAyIjI/nkD5GD6o94NwymevfuLUtoW3/U19fX\nV7ZtyJAhzWqX9u3bh0mTJtl9HLV9pqf+QkJCkJqaCo1GgxkzZjT77kFjCWoLCwtx7tw5KdHmnj17\nEB0djezsbKhUKowYMUJKTvvee+855fWZNDdBrb0YTDlJSEgItFotMjIyEB8f7+rqOEV6ejpOnTqF\npKQkeHt7AwBKS0uRkZEBhUKB8PBwAAyiiFqKtbmYBoNBdktEp9NJ8z4artBQUlKCu3fvSn97e3vL\nRpqGDh0qfb69vb0xdOhQAMDVq1dx9epVfp47OD8/P6xevRp5eXlITk5GUFAQZsyYYXfwbUpQO2rU\nKLMEtX/84x9x9OhRvPjiiwCsJ6fdunUrhg8fjuPHj6NHjx5OeX3NSVDbHJyA7kRGoxErV65EYGBg\nu8+ZZDAYkJiYiGXLlpllm/3hhx/w6aef4uc//zmOHj3KIIqonTAajbLpCHq9Xrp1X1VVJdt2/vx5\nab6JiUqlkn5vOPpVfzSiYcDWUeh0OigUCrNRwY7CaDQiPT0dx48fR2RkpNRhtsXOnTvx5ptvYvr0\n6dixY4dULoTAhAkToFAoWj7XkwXjx4+HQqGAVqs1y6vmzAnoHJlyIqVSiUWLFiExMRG+vr7tOrjQ\narUYM2aMxbT9TzzxBHr37o3s7Gy88cYb7fp1ErkTpVIpC3oau6Uza9YsrFq1SgqK6o92mUbGTF9C\nxcXF0naj0YgLFy7IzuXn5yeNfimVSlmg5efnJxsFqR+wtRVarRZpaWnS8iY//PADgoODO8xdCBPT\nWn86nQ7p6enQarWIiYmxKTBuKkHt888/j8zMzFYdaDAlqD18+LAskLp8+TJOnDjh1GsxmHIyX19f\nrFixAklJSRg6dGi77Z0ZjUb4+PhY3d6jRw+EhoYykCLqgNRqNXr37i1rv5o7l6bh6FdJSYkUhOn1\nelmQ5op5YU3R6XTYtWsXEhMTpc5lVVUVPvvsM6xduxYJCQlOv6armdb6y8zMRFJSEkJDQxEREdHo\nv6+1YAp4mO08JycHcXFx6NOnD4KDg1us7ib5+fn44IMPkJubi6eeekq2zZSl3Zl4m6+FqNVqZGRk\nYPny5e0yoMrMzMS9e/fw6quvWty+YcMG/OxnP2OqA6IOaOXKlQgODnZpZ8nSvLD62+ovJ9LceWGA\necDW0Pz58zF//nyzUfqqqiqsWLECixYt6tDtoMFgQFpaGgwGA2JiYqxO5O7Xrx+uXr2KK1euoG/f\nvhb3uXDhAhYuXIisrKyWrDKAh3O1NmzYIFssuSUxmGpB7TmgKiwsxNatW/Gb3/zG4vYVK1ZgxYoV\nHXLeAJE7M82XTElJaZdPJtszL6xhUAbIbzN2794dZ86cwX/9139ZvFZOTg66du3a7ufI2kKr1SI9\nPR1Dhw7FjBkz2PY3wGCqhanVamzevBmzZ89uF7fEtFot8vLyUFJSAqVSiUmTJuG5556T7ZOVlYWe\nPXu6RQNC5G4yMzNhMBjcMuN5w9EvnU6H0tJSLF682OL+RUVF2L59O6qrqxudF9bwlmRbnBdmC6PR\niKysLGg0GrsnqHd0DKZagV6vR1JSEiZPnozIyEhXV8eMwWCARqOBWq0G8HB41DSUu3btWnTp0kW6\nx/zdd99h6NChiImJaZe9ViJq3Pz58xu9neNOdDodtm/fjsTERIvbc3JyUFtbi9DQULN5YSZ6vR5V\nVVXS3/bMC2t4S7Kl5oXZS6/XIzU1FQqFwqHcVB0Jg6lWYjAYkJycDKVSiYSEhDbxgTBlNddoNHjx\nxRcREhJi8UOh1Wqh0+ng7e2NkSNH8oND1EHp9XokJycjJSXF1VVpM+Li4mS59urbsGEDnnnmmWaP\n0LTUvDBAHpS1VCqHvLw8ZGZm2jRBvaNjMNXK0tPTodFoEBMT0ypPNDSk1+ulUShvb2+Eh4cjODjY\nrT8ERPRQamqqlBmbHsrLy0NZWRkiIiJk5Tk5OTh16hRWr17d6nVqbF4YIB/9ampeWMNbkPVTVdiS\nL8xoNCItLQ3FxcWIj4932xFNBlMuoNPpkJaWBqVS2SpDpKYAypQJNjAwECEhIe1uUjwRtSxTbilO\nLpZLT0/H+fPn4e/vDwA4e/Ys7t27h3fffbfdtaOW8oWZ1L8laW++sAcPHkCtVqNnz5547bXXzOba\ndnQMplzINEQ6dOhQhISEOG2kymg0ori4WLo9BzzsiQQFBbltr4GIGqdWq5GXl+eSkZb2QKfTQafT\noby8HAEBAQgMDHSrEX1L+cJMTKNftbW1uHHjhtlIWHucF2YvBlMuZjQapSforl27hpCQEKhUKvj5\n+dncOzQYDCguLoZer0dJSQmKi4sxZMgQBAYGIigoqN31nIio9SUnJ0uj1kSO0Ov1SE9PR1VVFX76\n05+iZ8+e0rb6I2MNU1W0xXlhtmIw1Ybo9XoUFBSguLgYxcXFsjdH/TeNaWjW9Eb08vLC0KFDoVKp\npP+2x8ieiFzDYDBgwYIF+Oyzz9h2kNM4MkG9JeeFtcRT9Qym2jDTm6d+JG9iitAbvkmIiOyVl5cH\nnU6HRYsWuboq1MEYDAakp6e36gT1puaFtcT7nMEUAQCOHTuGpKQkHD58GBUVFbJtznqL2LuUQFRU\nFNavX29xOYCGq3/zbUzUfMwtRS2to2dQ7+TqCpDrHTp0CC+88AIGDRokm1S4d+9elJWVyfZNTU3F\nwIEDoVAorP706dPH7Br5+fkICwvDkiVLbK7X4sWLERoaim+++cZsW1lZGfbu3Wv7iyQii0wj4O4U\nSDVssyw5duwYpkyZAh8fH3h6emLMmDG4cOGCXbeIoqKiUFpa6lA9OorAwECsWrUKvXv3xrJly5CX\nl+fqKjmXILdWW1sr/P39hZeXl7h9+7YQQggAorG3RkVFhQAg+vbtK5XV1NSIqqoq8f3334uwsDDZ\n/t9++63w9vYW2dnZdtcvKytLKJVKUVhYaHF7U3UlosZt2rRJbN++3dXVcMimTZvEgAEDpPbA0o+v\nr6+0v6ls7969oqyszOx8Bw8eFJ07dxazZ88Wt27dEkeOHBEvvPCCePzxx8WRI0dsrteRI0eEv7+/\n0Gg0FreXlZWJvXv3ul07VlJSIpYuXSoSExNFSUmJq6vjFO7zf48sOnDggAAgwsPDpbKmPtj5+fkC\ngHj55ZfNtl28eFHMmzdP+ruurk6MGTNGBAUFNbuOgYGBYuzYsaKurs5sm7s1QkTOVFlZKeLi4kR5\nebmrq+Iwezp5jbUbljqYLdkhbKo+HVlubq6YOXOmSE9PF5WVla6ujkN4m8/NHThwAAAQGhpq8zFF\nRUUAgDFjxphtGzRokGwpitzcXBQVFSE6OrrZdYyOjkZhYSH27dvX7HMQkTmtVmtXGpa2zNQujR49\nWip75JFH4OXlBaVSiVGjRtl0noMHD+Lf//43QkND0b17dwghEBsbi4CAAEybNs3uekVEREClUiE2\nNpZzOxsIDw/H6tWrUVJSgmXLlll82Kq9YDDlprRaLYKCgrBmzRoAwNKlS63Od2qosLAQgLzRAoDX\nXnvNbN/s7GwAwNixY2XlJ0+eRP/+/aFQKLBjxw4AwO7du6FQKKRFlU1M1/niiy9seWlEZKN9+/a5\nZFmrlmBPJ68xDTuY7BC2LF9fXyxfvhwxMTFITk7G2rVrZSkQ2gsGU24qMDAQBQUFUuBy8eJFCCFQ\nXl7e5LGWeoBlZWWora012/f48eMAgICAAFn56NGjcfLkSXTu3FkKuDw9PaFQKPDRRx/J9jXl2DKd\ni4gcp9frUV5e3mGSdNrTybPEWgczKioKgHmHELC9U8gOYdMCAwORkpKC3r17Y8GCBe1ugjqDKTd2\n48YNlJaWwsfHx2L6AUtqampw+vRpAA8bF9NTKAMHDrTY2JjWdnr00UfNtvn6+mLy5MnYv38/zp8/\nj/j4eGzZsgUzZsyQ7dejRw8AwPnz5+16fURkXW5ubocJpAD7OnmWWOtg/uhHPwJg3iE0XcuWTiE7\nhLYxLbKdkJAAtVqNlStXyhJ3tmUMptxYY8Pi1pw5cwbV1dUICAiAePgAA4QQmDNnDsaNG2e2/507\ndwAAXl5eFs8XGxuLqqoqvPTSS0hJSUFsbKzZPqasuaZzEZFjDAYDvvnmG4SHh7u6Kk5hbyfPGksd\nzMY6hIBtnUJ2CO0zcuRIrF69GiqVCsuWLUNWVlabv/XHYMqNNSeYMg2lN2ygPDw8LDZaptwp1nKo\nmHp9KpUKU6dObfTaHT0PC1FrUavVePHFFzvExHPA/k6eNZbaxKY6hEDTnUJ2CJsnKioKGzZsgE6n\na/MT1BlMubHmBFPWjklJScHw4cPN9u/WrRsA4P79+2bbTpw4gRkzZmDq1Kk4cOCALOV/fffu3QMA\ndO/e3eZ6EpFlRqMR+/btk+YCdQT2dvKssdS+NdUhBGzvFLJDaD/TBPXIyMg2PUGdwZQbc2YwBTwc\nav/JT34iKzOt+H39+nVZ+ddff42IiAhkZGRg7ty5ePDgAXbu3GnxmqYFLIcNG2ZzPYnIsqysLKhU\nqg4zKgXY38mz5zyNdQgB2zqF7BA6LiQkBCkpKfDy8mqTE9QZTLmpmpoa6HQ6eHh4SJMjbWGpsRFC\n4O7du/jTn/4kBU8mpiUqzp07J5Xt2bMHkydPxsaNG6FSqRAWFgYfHx9s27bN4jVNS9y403IXRC3B\nYDBArVYjJibG1VVxKns7efacx1qHELC9U8gOoXMolUrMmTOnTU5QZzDlpr7//ntUV1fjySefRJcu\nXZrcf9GiRRgwYACuXr0KAOjXr580ybNTp07w9vZGfHw8wsLCZMeZ1rE6evSoVDZr1izcv39fymnl\n4eGBqVOnoqioCMOHD8etW7dk5ygoKAAAvP76681/wUQEX19frFixokONSgH2d/IssdbBtNQhBOzr\nFLJD6Fz1J6gnJSUhIyPD9bf+XJN4nVxtx44dAoCYPn262TY4cWmDuro6MW7cODF+/Phmn8N0PJeT\nIaL6EhISRP/+/Rtdkw+AyMjIkB1nqd347rvvBAAxatQoWXleXp4AID755BNZua+vrwAgjh07JpW9\n9dZbAoDw9/cXN2/elMo//vhjAUDs27fP4utgO9Z85eXlIikpScybN08UFBS4rB78v+emlixZIgCI\njz/+2Gyb6YNtbRFQexUVFQmlUik+//xzu4/961//Krp16yZOnjwpK3fXBUKJyHGW2g1rHcyW7hBa\nqw/Zp6CgQMybN08kJSW5ZK1J3uZzU7ZMPp8yZYrZ0i7NMXr0aOTk5CAxMREajcbm4/Lz8/HBBx8g\nNzcXTz31lGzbgAEDMGXKFIfrRkTuKycnB5cvXwZgvU1UKBTYtm0bzpw5g8zMTLuvsWvXLpw9exbb\ntm0ze5rv8uXLyMnJaWbtqb7AwECsWrVKlpuqNSmE4MqL7qhfv364evUqrly5gr59+7bKNS9cuICF\nCxfa/CaPjIzEhg0bbM7OTkRki4ZBjRACEydOxFdffYX9+/fjlVdeMTtGrVYjLi4OW7ZssXk9w/z8\nfMyaNQtbt27Fiy++aFM9yHF6vR7p6ekwGAyIj4+36yGr5mIwRUREbs+WDiY7hO2LWq1Geno6AgIC\nMGPGjBZ98ILBFBEREXVIRqMRWVlZ0Gg0CA8Px6RJk6SM9M7EYIqIiIg6NIPBgLS0NBgMBkRFRdl8\nq9ZWDKaIiIjILWi1WqSnp8PX1xfx8fFOu/XHYIqIiIjcilqtRkhIiNPOx2CKiIiIyAHMM0VERETk\nAAZTRERERA5gMEVERETkAAZTRERERA5gMEVERETkAAZTRERERA5gMEVERETkAAZTRERERA5gMEVE\nRETkAAZTRERERA5gMEVERETkAAZTRERERA5gMEVERETkAAZTRERERA5gMEVERETkAAZTRERERA5g\nMEVERETkAAZTRERERA5gMEVERETkAAZTRERERA5gMEVERETkAAZTRERERA5gMEVERETkAAZTRERE\nRA5gMEVERETkAAZTRERERA5gMEVERETkAAZTRERERA5gMEVERETkAAZTRERERA5gMEVERETkAAZT\nRERERA5gMEVERETkgP8PvkQF9fQ2AAAAAklEQVQb0r0QxdgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import Image\n",
    "\n",
    "Image('./jensens-inequality-1D.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\mathrm{log A} \n",
    "&= \\sum_{i=1}^{m} \\mathrm{log} \\bigg (\\sum_{z^{(i)}} p(x^{(i)}, z^{(i)} | \\theta) \\bigg) + \\mathrm{log}(p(\\theta)) \\\\\n",
    "&= \\sum_{i=1}^{m} \\mathrm{log} \\bigg (\\sum_{z^{(i)}} Q_i(z^{(i)}) \\frac{p(x^{(i)}, z^{(i)}| \\theta)}{Q_i(z^{(i)})} \\bigg) + \\mathrm{log}(p(\\theta)) \\\\\n",
    "&\\geq \\sum_{i=1}^{m} \\bigg (\\sum_{z^{(i)}} Q_i(z^{(i)}) \\mathrm{log} \\frac{p(x^{(i)}, z^{(i)}| \\theta)}{Q_i(z^{(i)})} \\bigg) + \\mathrm{log}(p(\\theta)) \\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $Q_i$ is a distribution over $z$ for data point $i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the E-step, we want to tighten the lower bound (i.e. convert $\\ge$ into $=$), which is easily achievable by setting \n",
    "$\\frac{p(x^{(i)}, z^{(i)}| \\theta)}{Q_i(z^{(i)})}$ a constant, i.e. let\n",
    "\n",
    "$$\\frac{p(x^{(i)}, z^{(i)}| \\theta)}{Q_i(z^{(i)})} = c$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, since $Q_i$ is a probability distribution, so it should sum up to one over $z$. So,\n",
    "\n",
    "\\begin{align*}\n",
    "\\sum_{z^{(i)}} Q_{i}(z^{(i)}) = \\sum_{z^{(i)}} \\frac{p(x^{(i)}, z^{(i)}| \\theta)}{c} = 1\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then,\n",
    "\n",
    "$$ c = \\frac{1}{\\sum_{z^{(i)}} p(x^{(i)}, z^{(i)}| \\theta)} = \\frac{1}{p(x^{(i)}| \\theta)} $$\n",
    "\n",
    "So,\n",
    "\n",
    "\\begin{align*}\n",
    "Q_i(z^{(i)})\n",
    "&= \\frac{p(x^{(i)}, z^{(i)}| \\theta)}{p(x^{(i)}| \\theta)}\n",
    "\\\\\n",
    "&= p(z^{(i)} | x^{(i)}, \\theta)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EM algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, we could formalize the EM algorithm for MAP as below,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the **E-step**, for each data point $i$, set\n",
    "\n",
    "$$Q_i(z^{(i)}) = p(z^{(i)} | x^{(i)}, \\theta)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the **M-step**, parameterize $\\theta$ to maximize $\\mathrm{log A}$,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\theta &= \\arg\\!\\max_{\\theta} (\\mathrm{log A}) \\\\\n",
    "&= \\arg\\!\\max_{\\theta} \\bigg[ \\sum_{i=1}^{m} \\bigg (\\sum_{z^{(i)}} Q_i(z^{(i)}) \\mathrm{log} \\frac{p(x^{(i)}, z^{(i)}| \\theta)}{Q_i(z^{(i)})} \\bigg) + \\mathrm{log}(p(\\theta)) \\bigg] \\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that RHS is a particular form of $\\mathrm{log A}$ after applying Jensen's inequality, I would like to label it as $Z$ for convenience later. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspecting this algorithm, what it's doing is basically to \n",
    "\n",
    "1. estabilish a lower bound of $\\mathrm{log A}$ in the E-step, and then\n",
    "1. tune parameters ($\\theta$) in order to maximize the lower bound."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Detour 2**: about the name Expectation-Maximization:\n",
    "\n",
    "Expectation probably refers to the quantity:\n",
    "\n",
    "$$\\sum_{z^{(i)}} Q_i(z^{(i)}) \\mathrm{log} \\frac{p(x^{(i)}, z^{(i)}| \\theta)}{Q_i(z^{(i)})}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm is indeed about maximizing the sum of expectations over all data points, but I thought the name might be too generic and not illustrative of the process of lower bound estabilishing. Anyway, I interpret the name as following:\n",
    "\n",
    "* E-step is the step of setting an expectation\n",
    "* M-step is the step of maximizing the expectation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next, we will prove each EM iteration is monotonically improving $\\mathrm{log A}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose after an EM iteration at step $t$,\n",
    "\n",
    "* $Q_i^{(t)} = p(z^{(i)} | x^{(i)}, \\theta^{(t-1)})$\n",
    "* $\\theta = \\theta^{(t)}$\n",
    "\n",
    "so"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\mathrm{logA}(\\theta^{(t)}) \n",
    "&= \\sum_{i=1}^{m} \\mathrm{log} \\bigg (\\sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) \\frac{p(x^{(i)}, z^{(i)}| \\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} \\bigg) + \\mathrm{log}(p(\\theta^{(t)}))\n",
    "\\\\\n",
    "&\\ge \\sum_{i=1}^{m} \\bigg (\\sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) \\mathrm{log} \\frac{p(x^{(i)}, z^{(i)}| \\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} \\bigg) + \\mathrm{log}(p(\\theta^{(t)}))\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we start a new EM iteration at step $t+1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the E-step, we set a new lower bound for each $i$ with\n",
    "\n",
    "$$Q_i^{(t+1)}(z^{(i)}) = p(z^{(i)} | x^{(i)}, \\theta^{(t)})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then replacing $Q_i^{(t)}$ with $Q_i^{(t+1)}$ in the above equation,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\mathrm{logA}(\\theta^{(t)}) \n",
    "&= \\sum_{i=1}^{m} \\mathrm{log} \\bigg (\\sum_{z^{(i)}} Q_i^{(t+1)}(z^{(i)}) \\frac{p(x^{(i)}, z^{(i)}| \\theta^{(t)})}{Q_i^{(t+1)}(z^{(i)})} \\bigg) + \\mathrm{log}(p(\\theta^{(t)}))\n",
    "\\\\\n",
    "&= \\sum_{i=1}^{m} \\bigg (\\sum_{z^{(i)}} Q_i^{(t+1)}(z^{(i)}) \\mathrm{log} \\frac{p(x^{(i)}, z^{(i)}| \\theta^{(t)})}{Q_i^{(t+1)}(z^{(i)})} \\bigg) + \\mathrm{log}(p(\\theta^{(t)}))\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the first equlity applies with whatever distribution $Q_i$, and the second inequality becomes equality become we reached the lower bound with the newly set $Q_i = Q_i^{(t+1)}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then in the M-step, we maximize ${\\mathrm{log A}}$ to obtain $\\theta = \\theta^{(t+1)}$, so"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\mathrm{log A}(\\theta^{(t+1)})\n",
    "&= \\sum_{i=1}^{m} \\mathrm{log} \\bigg (\\sum_{z^{(i)}} Q_i^{(t+1)}(z^{(i)}) \\frac{p(x^{(i)}, z^{(i)}| \\theta^{(t+1)})}{Q_i^{(t+1)}(z^{(i)})} \\bigg) + \\mathrm{log}(p(\\theta^{(t+1)}))\n",
    "\\\\\n",
    "&\\ge \\sum_{i=1}^{m} \\bigg (\\sum_{z^{(i)}} Q_i^{(t+1)}(z^{(i)}) \\mathrm{log} \\frac{p(x^{(i)}, z^{(i)}| \\theta^{(t+1)})}{Q_i^{(t+1)}(z^{(i)})} \\bigg) + \\mathrm{log}(p(\\theta^{(t+1)})) \\\\\n",
    "&\\ge \\sum_{i=1}^{m} \\bigg (\\sum_{z^{(i)}} Q_i^{(t+1)}(z^{(i)}) \\mathrm{log} \\frac{p(x^{(i)}, z^{(i)}| \\theta^{(t)})}{Q_i^{(t+1)}(z^{(i)})} \\bigg) + \\mathrm{log}(p(\\theta^{(t)})) \\\\\n",
    "&= \\mathrm{log A}(\\theta^{(t)})\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here,\n",
    "\n",
    "* the 1st equality is just the definition of $\\mathrm{log A}$\n",
    "* the 2nd inequality is again an application of Jensen's inequality\n",
    "* the 3nd inequality is a result of the maximization process, in which $\\theta^{(t)}$ becomes $\\theta^{(t+1)}$. Note, $\\mathrm{log}p(x, z|\\theta)$ and $\\mathrm{log}p(\\theta)$ are supposed to be concave in $\\theta$ and tractable as assumed in the lecture notes. So the maximization process is tractable\n",
    "* the 4th equality is just a restatement of the last equation after $Q_i^{(t)}$ => $Q_i^{(t+1)}$ replacement"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, we proved, EM monotonically increases $\\mathrm{log A}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: The derivation is actually very analogous to that in the lecture notes on EM: http://cs229.stanford.edu/notes/cs229-notes8.pdf except for the additional standalone $\\log{p(\\theta)}$ in the log likelihood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
